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The  power  series  method  used  by  the  author  to  generate  highly  accurate  finite 
difference  schemes  for  ordinary  differential  equations  in  II]  and  for  the  heat  equa- 
tion in  [2]  is  here  applied  to  the  wave  equation.  The  emalysis  runs  parallel  to  [2] 
and  involves  semi -discrete  approximations  in  t and  in  x before  the  totally  dis- 
crete scheme  is  derived.  The  results  differ  from  [2]  in  that  cm  arbitrarily  accurate 
difference  scheme  is  found  for  the  wave  equation  that  is  stable  and  consistent  with 
the  differential  equation.  No  such  scheme  exists  for  the  heat  equation.  The  step 
sizes  in  x and  t must  be  equal  for  this  difference  scheme.  Other  difference 
schemes  that  do  not  restrict  the  step  sizes  are  staible  only  when  the  order  of 
accuracy  in  x is  less  than  5.  The  lowest  order  scheme  is  shown  to  coincide  with 
Keller's  Box  Scheme  [3] . 

AMS  (MOS)  Subject  Classifications;  35A40,  65M05,  65M10. 

Key  Words:  Partial  differential  equation,  wave  equaltion,  power  series,  difference 
scheme,  high  accuracy. 
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SIGNIFICANCE  AND  EXPLANATION 


This  paper,  a sequel  to  TSR  #1923  and  TSR  #1924  where  the  author  introduced  the 
"Power  Series  Method"  and  applied  it  to  the  heat  equation,  applies  the  method  to  the 
wave  equation.  We  derive  difference  equations  that  approximate  the  differential  equa- 
tion by  substituting  power  series  into  it,  and  let  these  series  apply  in  a grid  of 
cells  that  cover  the  domain.  After  the  coefficients  in  the  series  are  evaluated  the 
series  are  truncated  at  some  order  of  accuracy.  As  much  of  the  detail  of  the  method 
is  described  in  the  above  quoted  TSR's,  this  paper  omits  such  detail  and  moves  rapidly 
through  the  analysis.  The  results  are  surprisingly  different  for  the  wave  equation 
compared  to  those  of  the  heat  equation.  Here  we  obtain  a set  of  difference  schemes 
of  increasing  accuracy  that  are  stable  to  all  orders  while,  for  the  heat  equation, 
only  schemes  accurate  to  order  15  in  x are  stable. 


POWER  SERIES  METHODS  III  - THE  WAVE  EQUATION 
Robert  D.  Small 

Introduction 

In  previous  papers  of  this  series,  [11,  [2],  which  we  refer  to  as  I and  II,  the  power 
series  method  has  been  used  to  derive  finite  difference  schemes  for  a general  first  order  or- 
dinary differential  equation  and  for  the  heat  equation.  This  paper  continues  the  study  of  the 
method  by  application  to  the  wave  equation.  Since  we  shall  depend  heavily  on  the  previous  work, 
a I or  II  preceding  equation  or  figure  numbers  will  refer  to  equations  or  figures  in  these 
papers.  In  addition,  since  the  treatment  of  the  wave  equation  follows  closely  that  of  the  heat 
equation,  a good  deal  of  detail  in  the  procedures  is  omitted  in  this  presentation.  As  in  II, 
we  discretize  the  wave  equation  first  in  t,  then  in  x and  finally  combine  the  two  discret- 
izations, maintaining  the  stability  conditions  that  have  accumulated  along  the  way.  It  turns 
out  that  one  can  obtain  a stable,  arbitrarily  accurate  difference  approximation  for  the  wave 
equation  whereas  this  is  impossible  for  the  heat  equation.  As  in  the  case  of  the  heat  equation 
the  root  locus  diagram  of  equation  (1-16)  illustrated  in  Figure  I-l  plays  a fundamental  role 
in  the  stability  analysis. 
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1.  Formulation  of  the  Problem 

We  now  prepare  to  derive  difference  schemes  for  the  wave  equation 


(1) 


■.2  .2 

3 u 3 u 


= 0 


3x 


3t 


The  domain  is  taken  identical  to  that  used  in  II  and  is  defined  by  0 ^ x ^ 2Lh,  0 ^ t < “>.  It 
is  partitioned  into  a grid  of  rectangular  cells  of  size  2h  by  2k  in  the  x and  t direct- 
ions respectively  and  each  is  identified  by  an  ordered  pair  of  integers  (i. j) . A local  origin 
is  located  in  each  cell  a distance  2Xk  from  the  bottom  of  the  cell  on  the  vertical  mid-line 
and  the  local  co-ordinates  are  denoted  (x^.tj , the  solution  to  (1)  then  being  denoted  within 
the  cell  by  u^^(x^,t^).  As  in  II,  since  the  preponderance  of  coefficients  requires  that  sub- 
scripts be  reserved  to  distinguish  them,  partial  derivatives  will  always  be  written  with  the 
symbol  3 and  the  superscript  (n)  will  denote  the  nth  ordinary  derivative  of  a function  with 
respect  to  its  argument.  Commas  are  used  to  separate  subscript  expressions  and  are  omitted 
when  no  ambiguity  arises. 

We  specify  that  initially  u is  U(x)  and  ^ is  V(x) . For  boundary  conditions,  u is 

o u 

g^(t)  on  the  left  side  of  the  domain  and  on  the  right.  Expressing  these  conditions  in 

local  co-ordinates , we  have 


(2) 


Uii(Xi,-2Xk) 

= U,(x.) 

3u . - 

3t' 

D 

= V,(x.) 

u,.(-h,t.) 

= 

u^.(h,t.) 

i = 1,2,- 


■,L 


j = 1,2,- 


3u 


For  interior  cell  boundaries  we  specify  that  u and  ^ are  continuous  across  horizontal 

boundaries  and  that  u and  ^ are  continuous  across  vertical  boundaries.  These  conditions 

3x 


1 


2.  Discretization  in  t 


Substitution  of  a power  series  in  t into  (1)  gives  the  local  solution 


Ui.(Xi,t  ) = I 


jj + y ^ 

'’-''I  -i  (2n+l)!  j 


j'  (2n)!  3 


Substitution  of  this  into  (2)  and  (3)  leads  to  the  following  infinite  system  of  ordinary  differ- 
ential equations: 


b‘fUi. 


V 3-1  1 , .zn+1  „ , . 

^ (2n+l)! 


,(2n+2),  , u(2n),  , 

(2n+l)!  (-2>k)  + I (-2Xk) 

n=0  n=0 


i = 1,2, •••,L 


(2n)  , (2n),„, 

a^j  (-h)  = g^.  (0) 


, (2n)  , . , (2n+l)  , 

(-h)  = gj^j  (0)  n = 0,1,  • • • ,“ 


a^f.h,  =g^f><0) 


j = 1,2,- 


b‘^"’  (h)  = (0) 

lij  ^ J 


(2n)  , , 

a,.  (x.) 


v,(2n),  , 

b. . (x. ) 


y ^ ^ (-2Xk)^*'  + y ^ j-  (-2Xk) 

(2n)!  ^ ^ (2n+l)!  ' 


aP"’,  (xj 


A. 

(2n)  ! 


, (2n) 


(2(1-X)k)  " t I (2(1-X)k) 


af^""2>(x  > 


b'^">(x.) 


'i  = 1,2,*'-,L 


j = 2,3,*'-,'> 


(2n+2)  , , 

” a.  . , (X.) 

y i- 

(2n+l)! 


n„'  v,<2n)  , , 

2n+l  ^ V ‘^i,j-l  ^i 


(2(1-X)k)^""'‘  + I 


(2(1-X)k)' 


As  is  usual  in  the  semi-discrete  formulation,  (5d)  indicates  that  the  partitioning  into 


cells  in  the  x-direction  is  unnecessary;  the  problem  takes  place  in  horizontal  strips.  Equa- 
tions ( 5c)  are  truncated  in  the  following  manner . 

”2  ""l  ”2 

I - I = I ^ I 

n=0  n=0  n=0  n=sO 

N3  N3  N4 

I - I = I ^ I 

n=0  n=0  n=0  n=0 

This  system  is  tested  for  stability  by  the  substitution 

. lax. 

aij(Xi)  = B^e 

. lax. 

b.  . (x.)  = B^e  ^ , 

where  i = 

The  result  is  a system  of  homogeneous  equations  for  which  the  condition  for  non-trivial  solu- 


tions is 
N 


r (-1)"[X^"  - B~^(l-X)^"  (2ak)^”  r (-1)”[X^"  - B~  ^ ( 1-X ) ( 2ak)  ^ 

i (2n)  ! ^ (2n)  ! 


n=0 

N 


n=0 


(2n)  ! 


I 

n=0 


(2n+l)  ! 


I 

n=0 


(2n+l) ! 


= 0 . 


The  choice  N,  = N.  and  N,  = N,  allows  scraare  roots  to  be  taken  and  also  suppresses  some  un- 
14  2 3 

N-1 

desirable  roots  for  6 . The  the  best  accuracy  occurs  for  = [N/2]  and  = [“^1  fo’^ 
some  order  of  accuracy  N.  Making  these  choices  we  solve  for  B obtaining 

r (-l)"(2(l-X)ak)  ^ f (-l)"(2(l-X)ak) 

^ (2n) ! ^ ^ 


B = 


n=0 


n=0 


(2n+l) I 


I 


^ (-l)"(2Xak) 


n=0 


(2n)  ! 


-i  I 


^ (-l)"(2Xak)^"'*'^ 


n=0 


(2n+l) t 


and  its  complex  conjugate.  We  can  obtain  iB|  = 1 by  choosing  X = 1/2  but  other  values  of 
X do  not  consistently  maintain  |b1  ^1  for  all  a or  k,  hence  we  tcike  X = 1/2.  With 


-5- 


these  choices  we  select  the  correct  boundary  conditions  from  (5b)  and  (5d)  to  suit  the  differ- 
ential equations  and  obtain  the  followinq  semi-discrete  scheme  that  is  still  continuous  in  x. 

f "l  a‘2n)u.)  , ^^2  bif>(x.)  , , " 

V tl  1 2n  _ r rl  1 . 2n+l  _ , > 


i = 


N (2n+2)  N_  . (2n)  , , 

2a.,  (x. ) _ , , 1 b. , (x. ) 

r il 1 . 2n+l  y il  1,  2n  „ , , 

■ nlo  ‘ io  (2n)  - = '^i<='i> 


a<^"N-h)  = n = 0,1, 

b{^"’(-h)  = n = 0,1,---,N^-1 

aLj"’(h)  =g2j"’<0’  n = 0,l,---,N2 

=g‘^"‘^^’(0)  n = 0,1, •••,»-! 

LJ  ZJ  i 


N ^(2n),  , N-  .(2n),  , 

1 a,,  (X.)  _ 2 b, ■ (X.)  _ 

„io  (2n).  (2n.l)( 


j = 1,2, 


N,  (2n)  , , N,  , (2n)  , , 

la...  (x.)  _ 2 b...,  (x.)  ^ . 

y + y ^’3~^ i_  .-1  = 12  •••L 

{2n)!  (2n+l)!  ^ 

nstn  nssfi 


y ^'3~^  ^ k^"  + y a 

n=0  n=0 

•^2  ^(2n) 

y ^3  1 ,^2n+l  r 1.J  1 2n 

nio  n=0  '2n)! 


= 2,3,--' 


N (2n+2)  N,  , (2n)  , , 

y i.j-1  i’  2n+l  y \,j-l‘^i^  2n 

nio  nio  J 

|^a["’(-h)  = af"^^^(h)  n = 0,1,  • • • ,2N2+1^  fi  = 2,3,---,L 

I b^"’ (-h)  = b["^^  .(h)  n = 0,1,---,2N^-iJ  I j = 1,2,-.-,» 

N,  (2n)  , , N.,  ,_(2n)  , , , . , 

. , ’ ,2n  ^ ^1  2n+l  f"  = 1-2, •••,L 

' ^r,  <2n)!  *^j  ^ (2n+l)  ! \ . , . 

■'  n=0  n=0  'j  = l,2,>--,oo 


i-l,j<‘'’ 

i"!  .<h) 


u..(x.,t)  = ) ^ t^"  + y 

^ ^ n=0  <2n)!  , 
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We  begin  with  the  power  series  in  x.  which  satisfies  (1) 


2in+l 


(2m+l)  ! 


Substitution  into  (2)  and  (3)  gives  the  following  system,  where  we  have  taken  X = 1/2 


(2m+l) 


-(2III+1) 

‘il 


■(2m+l) 

il 


(2m+l) 


2m+l 


(2m+l) 


(2m+l)  ! 


2m+l 


(2m+l)  ! 


(2m+l)  ! 


(2m+2) 


2m+l 


(2m+l)  ! 


•(2m+2) 

i-i.i 


2m+l 


We  truncate  (8d)  in  the  following  manner. 


M 


M, 


M. 


1 “2  1 2 

I - I = I ^ I 

m=0  m=0  in=0  m=0 


M 


M. 


3 4 3 4 

I ^ I - I ^ I 

m=0  m=0  m=0  m=0 


and  test  the  stability  of  this  new  system  by  the  substitution 


6t.  ^ 

— , 3 2iloih'' 

a.  . (t.)  = e -^e  a 
11  J 


6t. 


, 3 2iiah* 

b.  . (t.)  = e -"e  b 
13  3 


The  resulting  system  of  homogeneous  equations  has  non-trivial  solutions  if  the  following  equa- 
tion holds. 

M, 


r (Bh)^”'  r (Bh)^ 

^ /0m\ I L 


„_n  (2ro)i  f:  (2m)  i 
ni=0  m=0 


To  suppress  undesirable  roots  in  B we  take  M, 
to  be  taken.  Then  we  set  = [1 


"2 

y 

L 

m=0 

we  take 

= M4 

and  = 

(2m+l) J (2m+l) I 

in*0 


and  which  allows  square  roots 


order  of  accuracy  M.  The  factor 
2 


- c . With  these  choices  we  obtain 

M 

(9) 


reduces  to  - cot  ah  which  we  set  equal  to 


M, 


r (Bh)^*"  _ 3 V (Bh) 

^ (2m)! 


(2m+l) ! 


in=0  m=0 

This  is  equation  (1-16)  whose  root  loci  for  parameter  c were  originally  displayed  in  Figure 
I-l.  This  diagram  is  symmetric  about  the  real  and  imaginary  axes  and  is  later  reproduced  in 
Figure  1.  For  stability  we  must  have  Re  B £ 0 for  all  a and  hence  all  c.  From  Figure  1 
we  can  see  that  there  are  loci  that  lie  in  the  first  and  fourth  quadrants  for  M > 4 and  hence 
this  semi-discrete  scheme  is  stable  only  for  M £ 4.  In  this  way  the  x-discretization  suffers 
a severe  restriction  on  the  order  of  accuracy  of  difference  schemes  that  may  be  used. 
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Now  selecting  suitable  boundary  conditions  from  (8a)  and  (8c)  we  obtain  the  following 


semi-discrete  scheme,  continuous  in  t 


(2m+l) 


(10a) 


■(2m+l) 

il 


-(2m+l) 

'il 


(2m+l) 


2m+l 


(2m+l)  ! 


(2m+l)  ! 


2m+l 


(2m+l)  ! 


2m+l 


(2m+l)  ! 


(lOd) 


M-  -(2m+2) 


2m+l 


(2m+l)  ! 


M,  -(2m+2) 


(2m+l)  ! 


4.  Full  Discretization 


We  now  discretize  the  two  semi-discrete  schemes  (6)  and  (10)  with  respect  to  the  variables 
that  remain.  To  discretize  (6)  we  set 


oo  c 

r ijm  m 


00  d 

r lim  m 


. .(X.)  = I , b..(x.)  = ^ 

11  1 m!  1 11  1 "i!  1 


Substitution  into  (6)  leads  to  the  system 


Ml  N2  ^ j 

r il,m+2n  , 2n  r il,m+2n  . 2n+l  ..(m) 

I OnTT'  - I = “i  <°> 


N2  ^1  d 

r ’^il,m+2n+2  , 2n+l  . r il,m+2n  , 2n  ..(m) 
- I (2ntl).  ^_L  ’UnM"  = '^i 


m = 0,1, • • • ,” 


i = 1,2,**',L 


00  C 

r li,m+2n  . , .m  (2n) ^ , 

1 — (-h)  = gj^.  (0)  n = 0,1, '••,N2 

m=0  ”*  ■ ^ 

00  (J 

r lj,m+2n  , ,.m  (2n+l)  , „ 1 

> — (-h)  = g,  . (0)  n = 0,1,”*,N  -1 

m!  ll  1 


"Li,m+2n  ,m  (2n)  „ , 

— h = goj  (0)  n=0,l,*",N 

m!  ^2  ^ 


j = 1,2, 


Li ,m+2n  m _ (2n+l) 


(0)  n = 0,1,  • • • ,Nj^-l 


r ‘^ii  ,m-t2n  2n  _ r ij  ,m+2n  . 2n+l 
{2n)  ! (2n+l)  ! 


r ‘^i,i-l,m+2n  2n  r i,i-l,m+2n  . 2n+l 
(2n)!  (2n+l)! 


”2  *^1  d 

Y *^11  ,m+2n+2  , 2n+l  y ij  ,m+2n  2n  _ 

(2n+l)!  L (2n)! 


m = 0,1 , • • • ,<» 
i = 1,2,"**,L 
j = 2,3,  • • • ,<» 


1 d.  . 


Y i, 1-l,m+2nt2  . 2n+l  r i,i-l,m 
^ (2n+l)  ! (2n)! 


-l,m+2n  j^2n 


I (-h)"'  = I ^LZliidSiS.  h"*  n = 0, !,•••, 2N  +1 


[ (-h)”  = I -■'--4—'*'”  h*"  n = 0, !,•••, 2N  -1 

~ m!  m!  1 

m=0  111=0 


i = 2,3,*”,L 


•isl  0 •••00 

J “•  -L  f ^ / t 


V r n ,m+2n  m 2n  . ^ V n ,m-t-2n  _,in^2n+l  / 


. . . . . r r iT,m+<!n  m.  r r ii,m+2n  m.  2n+l  i 

lie)  u (X  t ) = I I x.t  + I ^ s 

in=0  n=0  m=0  n-0  vj  = 1,2,***, 


To  discretize  (10)  we  set 


l.At.)  = I ^t"  , b..(t.)  = I ^t"  . 

"3  3 „io  n!  j 13  1 ^to  3 


OO  d ^ 

r iin 


Substitution  into  (10)  gives  the  system 

nio  "=  " ' 


il,2m+n  = of^m+l)  m = 0,1,  • • • ,M  -1 

n ! 1 J- 


y fil^JnHjvn  (-k)"  = v!^"- (0)  m=0,l,'",M2 

n=0 

y ,.k)"  = vi'^^^’cO)  m = 0,1,---,M  -1 

n(  1 1 


,i  “ 1#2»***»Ij 


r ^li.2m+n  . 2rn  r li,2m+n  , 2m+l  (n) 

L “(k)T“ ^ ■ 4 (Ml  ^ = ^ij 


M,  " 

1 C.  ^ 


”2  d.. 


Li,2m+n  . 2m  . r Li,2m+n  , 2m+l  (ri),,>> 

ikrr-  ^ * 4 imt  ^ = ^23 


n = 0,1,*’*,“ 


j = 1,2,-*-, 


r *^13  ,m-t-n  , , ,n  r *^1,3-1, m+n  n 


I - I ^ 

n=0  n=0 


k m = O,!,-- *,2112+1 


y liiJSlii  (.k)"  = y fLiliJSiiik"  m = o,i,..*,2m,-i 

^ n I n ! 1 


f i = 1,2,  • • • ,L 

/ 

\ j = 2, 3,  • • • ,“ 


■ ’'4^ 


C ”2  d 

r ii,2in+n  , 2m  r ijf2m+n  , 2m+l 


”1  c ”2  d 

r i-l.i,2m+n  , 2m  . r i-l,i,2m+n  . 2m+l 

J.  (2m;i)  ! 

m=0  m=0 


”2  “ ”1  d 

V ^ii*2m+n+2  , 2m+l  r ii,2m+n  , 2m 
(2mfl)  - ^ L = 


”2  c ”1  d 

r i-l,i,2m+n+2  . 2m+l  . r i-l,1,2m+n  ,2m 

-(2m+l)! ^ J. 

m— 0 n^O  J 


■»  c ^2  “ d i = 1,2,  • • • ,L 

V V ii  ,2m+n  2m  n r r ij  ,2iiv+n  2mfl  n ) 


n = 0,1,  • • • ,“ 
i = 2,3,”-,L 
j = 1 ,2, • • • ,“ 


^ 3 „,to  n=0  ■ ‘2m)  !n!  i 3 n=0  ^ 3 I j = 1,2,  • • • ,-  . 

Systems  (11)  and  (12)  express  finite  difference  equations  tliat  are  written  in  different 
variables  and  truncated  at  different  stages  of  the  series.  To  select  common  variables  we  com- 
pare (lie)  and  (12e)  and  note  the  correspondence 

'^ij,2n  '^ij,2n  '^ij,2n  '^ij,2n+l 

’^ij,2n+l  '^ij,2n  *^ij,2n+l  '^ij,2n+l 

If  (12)  is  now  expressed  in  the  variables  of  (11)  , we  can  choose  all  the  finite  limits  of  summa- 
tion from  this  new  system  and  from  (11)  and  obtain  the  fully  discrete  scheme.  That  final  sys- 


1 2 d 
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”2  c “id 

V il,nH-2n+2  . 2n+l  r il,m+2n  . 2n  ■ ..(m) 
■n=0  n=0  <2")  ‘ ^ " 


m = 0, 1, • ' • ,M-1 


i = 1,2,'”,L 
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m ^ (2n) 

m!  Ij 


n » 0,1, 


(-h)™  = (0)  n = 0,1, -.-.N, 


m! 


[ Q 

Li  ,iiH-2n 

0 


= gjj"’ (0) 


n = 0,1, •• • ,N, 


i ‘^Lj  ,in-»-2n  j^iti 

:o  ™‘ 


= g<f^'’(0) 


n = 0,1, • • • ,N 


N 


(2n)  ! 
N 


ij  ,m+2n  2n  r ,in4-2n  v2n+l 

(2nn  (2n+l)! 

n=0 


y ‘^i,i-l,nn-2n  ^2n  y '^i,i-l,m^-2n  ^2n+l 
L (2n)l  ^ /2n+ni 


n=0 


n-0 


(2n+l) I 


1 d. 


ij  ,iiH2n-t-2  . 2n+l  y i j ,m+2n 

/2n4.i » 1 L ^2n,  • 


2n 


(2n+l) ! 
N 


n=0 


r ‘^i,i-l,m-»-2n-»-2  2n+l  r ‘^i,i-l,iiH2n  2n 

(2n+l)!  ^ (2n)! 

n=0  n*0 


I . o,l.  • • • ,2N,+1 

m!  ^ ml  2 


M c. 


m*0 


m! 


i * 


fiisSiiii  (.h)"*  = I h"  „ = o,i,...,2n,-i 

ml  _ m!  1 

m=0 


3 = 


M 1 C.  . 


M ”2  d. 


‘ ) ..0  nSo  ‘ J io  „io  ‘ J 


2n+l 


c.  . 


ii,M+N  M.  N 

— x.t.  (N  even) 


M!N!  i j 


d.  . 


ii,M+N-l  M.  N 
•‘J,-. x.t.  (N  odd) 


i = 1,2,""*,L 


j = 1,2,  • • • ,“■ 


The  final  term  in  (13e)  arises  because  the  highest  coefficients  computed  in  (13a)  to  (13d)  are 

c. . , and  d. . „ The  double  sum  in  (13e)  contains  one  term  that  is  not  confuted 

1],M+2N2+1  1],M+2Nj^-1 

and  this  must  be  deleted. 

We  now  test  system  (13)  for  stability.  To  examine  one  Fourier  component  of  the  error  we 

set 

2iiah- 

c.  . = e e c 

ijn  n 

. -j  2iiaha 

d.  . = Be  d 

i3n  n 


Substitution  of  this  into  (13c)  and  (13d)  gives  the  following  set  of  homogeneous  equations. 

r 


(14) 


} (1-B'S  2n. 

n=0 


^2  -1 
-I 

m+2n  (2n+l) ! 

n=0 


m+2n 


= 0 


^2  -1 


N, 


\ m = 0, 1, • • • ,M-1 


n=0 

M 

I 

tn=0 

M 


(2n+l) 


+ y i,2n> 

■m+2n+2  (2n)  ! ' 

n=o 


m+2n 


= 0 


, , , . m -2iah. 

y ((.-.IJ.  Lh"’a 

m!  m 


n = 0,1,**-.2N2+1 


*'*  / / T \ ^ ”2iQh.  * 

r ((-1)  -e ^ h d . =0  n = 0, !,•••, 2N  -1 

ml  m+n  1 

m=0 


The  matrix  of  coefficients  of  (14)  resembles  the  Sylvester  matrix  encountered  in  II  except 
that  this  one  has  four  sets  of  staggered  rows  instead  of  two.  We  note  that  if  its  determinant 
is  set  equal  to  zero,  we  obtain  an  algebraic  equation  of  degree  2M  in  B and  thus  we  see)r 
2M  solutions  for  B for  each  fixed  o.  The  substitution  that  collapses  the  equations  is 

c - z",  n = 0, !,•••, M+2N,+1 

n * 

d - ± z"*^,  n = 0,l,”*,M+2N,-l  . 

n 1 


Then  the  two  linearly  independent  equations  that  survive  are 
‘ n=0  ' 1-B-’  n=0 


and 
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2m+l 


(16) 


1 

I 


(2in)  ! 


2 

I 

in=0 


(hz) 
(2m+l) ! 


where  c = cot  ah.  Equation  (16)  is  of  degree  M in  z and  for  fixed  c gives  M solutions. 
Then  equation  (15),  being  linear  in  6 but  double  valued,  furnishes  the  full  2M  solutions  for 

B. 

For  non-trivial  solutions  for  the  errors  we  must  have  that  equations  (15)  and  (16)  hold 

1+g-l 

simultaneously.  When  |6|  =1  the  factor  r is  purely  imaginary  so  that  (15)  is  of  the 

1-6' 

form  (1-16) . The  root  loci  of  this  equation,  displayed  in  Figure  I-l,  are  reproduced  here  in 
Figure  1 for  the  variable  kz.  These  curves  divide  the  complex  plane  into  regions  where  solu- 
tions of  (15)  fall  for  |6|  < 1 and  for  |6|  >1.  Since  (15)  is  double-valued,  however,  any 
point  in  Figure  1 that  corresponds  to  a solution  for  |6j  <1  also  corresponds  to  another 

solution  where  |6|  > 1.  Thus  any  value  of  Icz  not  on  a curve  or  the  imaginary  axis  of  Figure 

1 corresponds  to  an  unstable  error.  Equation  (16)  is  also  of  the  form  (1-16)  and  thus  Figure  1 
represents  its  root  loci  for  the  variable  hz,  parameterized  by  c. 

For  stability  of  any  difference  scheme  (13)  of  order  (M,N)  we  must  have  non-trivial 
errors  only  for  |6|  £ 1.  To  investigate  whether  this  is  the  case  we  identify  the  roots  hz  of 

(16)  as  the  curves  corresponding  to  M and  the  imaginary  axis  of  Figure  1.  These  curves  are 

mapped  into  the  same  diagram  but  for  the  scale  adjusted  to  that  of  the  variable  kz  so  that  the 
value  of  6 in  (15)  can  be  ascertained.  If  the  mapped  curves  fall  only  upon  the  curves  corre- 
sponding to  N and  the  imaginary  axis  then  all  errors  occur  for  |6|  =1  and  the  scheme  is 
stable.  For  any  other  case,  some  errors  correspond  to  Ib]  > 1 and  instability  results. 

Clearly  scheme  (13)  is  stable  for  all  h,  )c,  N if  M £ 4 since  the  solutions  of  (16)  lie 
on  the  imaginary  cucis  and,  under  the  stretching  of  a scale  change,  remain  there.  This  result 
we  anticipated  from  section  3.  An  additional  set  of  stable  schemes  is  obtained  by  setting 
M = N and  h = )c  since  the  roots  of  (15)  for  |b1  =1  and  the  roots  of  (16)  then  coincide. 

In  this  way  we  obtain  stable  difference  schemes  for  (1)  of  arbitrarily  high  order. 
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5.  Truncation  Error  and  Consistency 

We  define  the  truncation  error  of  a computed  quantity  as  the  difference  between  the  com- 
puted and  the  exact  solutions  for  it.  Thus  the  truncation  errors  c..„,  d and  u. . are 

ijn  ijn  1] 

given  by 

(e) 

C.,  “C..  -C., 

ijn  ijn  ijn 

d.  , = d.  , - d.‘®^ 

I3n  ijn  ijn 


(e) 

u.  . = u. . - u, . 

11  ID  ID 

where  c.j  , d. , and  u. . are  the  computed  solutions  and  satisfy  (13)  for  finite  M and  N 

ijn  ijn  ij  ^ 

/a)  fc) 

while  c.  . , d.  . and  u.  , satisfy  (13)  for  M = “■,  N = “.  Substracting  the  two  forms  of  (13) 

ijn  i^n  ij  ■' 

we  obtain  the  following  equations  in  the  errors: 

r '*1  c ''2d  ^ 

p il,m->-2n  2n  r il,mt2n  , 2n+l  _ 
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) : I 


Me  ■»  c^®' 

p i j,m-f2n  _ y ii,m-t-2n 
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Li  ,m+2n  . m 

— — : h 
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If  the  determinant  of  the  coefficient  matrix  of  (17a)  to  (17d)  is  not  zero  then  the  solutions 
' " M+1  N+1 

for  the  errors  c. . and  d. . are  0(h  ,k  ) . If  this  is  so  then  it  follows  from  (17e) 
ijn  i^n 

- n+1  N+l  M N 

that  is  0(h  ,k  ,h  k ) . The  difference  scheme  (13)  is  consistent  with  the  different- 

ial equation  (1),  the  boundary  conditions  (2)  and  the  continuity  conditions  (3)  since  the  errors 
vanish  as  h and  k approach  zero,  for  all  M ^ 1,  N ^ 1,  provided  that  the  coefficient 
matrix  of  (17)  is  of  full  rank. 

This  last  condition  is  investigated  by  examining  the  homogeneous  system  formed  by  deleting 
the  terms  on  the  right  side  of  (17)  cind  the  errors  at  time  step  j-1.  Then  we  may  delete  the 


subscript  j obtaining  the  system 
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m! 


~ . . Ill  "s  . Ill 

M d.  ^ -h)  -d.  T h 

1 ,mtn 1-1  ,m4-n 


m=0 


m! 


=0  n = 0,1, • • • ,2N2+1 


0 n = 0,1,  • • • ,2Nj^-1 


\ i = 2,3, • • • ,L 


The  following  substitution,  motivated  by  the  treatment  of  (14),  turns  out  to  give  general 
solutions  for  (18): 


(19) 


i , 2n 

. 2n 
= A . z 

1 

i,2n+l 

„ 2n+l 
= B . z 

1 

‘i,2n 

^ 2n 
* C.z 

1 

'i,2n+l 

r.  2n+l 
= D.z 
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On  substitution  of  (19)  into  (18) , the  following  notation  becomes  convenient. 


«1  . 2m  »2 

liu  ' = I 


(hz) 
(2m+l)  ! 


N N 

^ 2 ,2n+l 

SO'-'"  I Sr  ' S0'->=  I SiTT 

n*0  n*0 


Then  sxibstitution  of  (19)  into  {18a)  results  in  a number  of  equations  that  differ  by  factors 
that  are  powers  of  z.  The  linearly  independent  equations  that  survive  are 


For  non-trivial  solutions  we  require 

(21)  K^()c,z)  = ± K^Oc.z) 


and  then  we  have 
(22) 


C.  = ± zA,  'I 

^ ^ I i 

D.  = ± zB.  j 

1 1 


1,2 


Placing  (19)  into  (18b)  emd  using  (20)  and  (22)  we  again  obtain  a large  set  of  equations  that 


For  (13)  to  be  consistent  with  (1),  (2)  and  (3)  we  must  not  have  non-trivial  solutions  for 


the  errors  in  (18),  These  errors  are  non-trivial  if  (21)  is  satisfied  and,  in  addition,  there 


are  non-trivial  solutions  to  (23) . To  test  for  consistency  for  a given  choice  of  M,  N,  L,  h, 
)c,  we  place  Ic  into  (21)  to  determine  z and  then  if  the  determinant  of  the  coefficient 
. matrix  of  (23)  is  non-zero,  (13)  is  consistent.  It  would  seem  unlUcely  that  a random  choice  of 

h and  )c  would  give  an  inconsistent  scheme.  We  found  in  section  4 that  to  obtain  a stable 
scheme  for  M ^ 5 we  were  forced  to  take  h = )c  and  M = N.  This  choice  happens  to  give  a 
consistent  scheme  since  (21)  leads  to  = ± and  then  the  rows  of  the  matrix  in  (23)  are 
mutually  orthogonal  which  implies  that  only  trivial  solutions  exist.  Thus  we  have  a stable, 
consistent,  and  arbitrarily  accurate  difference  approximation  for  (1) . 

As  in  II,  it  is  difficult  to  prove  that  substitution  (19)  yields  all  solutions  for  (18) 
since  one  does  not  have  a general  formula  for  the  determinant  of  the  matrix  in  (23)  . The  worlc- 
ing  of  several  examples,  however,  indicates  that  (19)  is  sufficiently  general. 

I 
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6.  Lowest  Order  Scheme  and  Concluding  Remarks 


We  have  indicated  in  II  that  the  lowest  order  scheme  generated  by  the  Power  Series  Method 

for  the  heat  equation  coincided  with  Keller's  Box  Scheme  [3]  for  the  heat  equation.  Making  the 

same  comparison  here  we  set  M = 1,  N = 1 in  (13)  and  eliminate  the  coefficients  c. . and 

ijn 

d^^^  in  favor  of  the  solution  values  obtain  the  following  difference  scheme,  arranged 

to  display  the  approximation  of  the  derivatives  of  (1) . 

(u.  . -2u.  , .+u.  , . + 2(u.  . ,-2u.  .+u.  . ,)  + (u.  , . ,-2u.  , .+u.  , . 

1-1, T+l  1-1,1  1-1, T-l 1,1+1  n 1,1-1 1+1, 1+1  1+1,1  1+1, 1-1) 

4(2k)^ 


(“■  , -.i-Zu.  +u  . ) + 2(u.  , .-2u.  .+u.  , .)  + (u.  , . ,-2u,  . +u.  , . ,, 

1-1, 1tl  1,1+1  1+1, 1+1 1-1,1  1,1  1+1,1 1-1, 1-1  1,1-1  1+1, 1-1) 

4(2h)^ 


Alternatively,  if  an  auxilliary  variable  v is  introduced  into  (1)  so  that  we  have 


3u  _ ^ ^ ^ 

3x  3t'  3t  3x 


then  the  Box  Scheme  (3)  used  on  the  derivatives  of  this  system,  with  the  resulting  varicibles 
V eliminated,  gives  the  above  difference  scheme.  This  demonstration  strengthens  our  specula- 
tion made  in  II  that  the  lowest  order  scheme  in  the  Power  Series  Method  applied  to  any  differ- 
ential equation  gives  the  same  result  as  Keller's  Box  Scheme  when  the  latter  can  be  applied. 

In  conclusion,  having  applied  the  Power  Series  Method  to  the  wave  equation,  arbitrarily 
accurate  difference  schemes  have  been  produced.  The  point  of  expansion  of  the  power  series 
within  each  cell  was  taken  as  the  center  of  the  cell  in  order  that  the  t-discretization  be 
stable.  The  heat  equation  at  this  stage  required  only  that  this  point  be  at  least  as  high  as 
the  center  point  and  was  later  chosen  to  be  the  center  point  to  make  the  total  discretization 
stable.  We  speculate  that  the  center  point  of  the  cell  should  always  be  the  point  of  expansion 
of  the  series.  The  stability  tests  of  sections  2 and  3 indicated  that  the  t-discretization  was 
unconditionally  stable  and  that  the  x-discretization  was  stable  only  for  M 4 . The  stability 
test  for  the  combine  discretization  corroborated  these  results  in  that  M £ 4 was  uncondition- 
ally stable  for  any  N.  Another  stable  set  of  difference  schemes  that  was  not  predicted  by  the 
semi-discrete  stability  tests  was  the  case  where  h = k and  M = N.  This  is  the  set  of 
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20.  AsSTRACT  (Contlnu*  on  rmvtoo  oldo  U nmeooomry  md  Idontlty  by  block  numbt) 

The  power  series  method  used  Jay  -the  autho^  to  generate  highly  accurate  finite 
difference  schemes  for  ordinary  liLfferen^i^l  equation^/in  [1]  and  for  the  heat 
equation  in  t2>-'is  h^e  applied  to  the  wave  equatibn.  The  analysis?  runs  paral- 
lel to  [21  anchJinvolves  semi-discrete  approximations  in'  t“  and  in  x before 
the  totally  discrete  scheme  is  derived.  The  results  differ  in  that  an 

arbitrarily  accurate  difference  scheme  is  found  for  the  wave  equation  that  is 
stable  and  consistent  with  the  differential  equation.  No  such  scheme  exists 
for  the  heat  equation.  The  step  sizes  in  x and  t must  be  equal  for  this 
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difference  scheme.  Other  difference  schemes  that  do  not  restrict  the  step  sizes  are 
stable  only  when  the  order  of  accuracy  in  x is  less  than  5.  The  lowest  order  scheme 
is  shown  to  coincide  with  Keller's  Box  Scheme^ [3]. 
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